% Shells for probability tables
prob2=[]; prob3=[];

%% Construct data for estimating VAR for y_t.
YY = data(opt.p+1:end,:); % y_t 
XX = lags(data,1:opt.p); % Matrix of regressors in VAR for y_t
XX = XX(opt.p+1:end,:); % Drop initial missing observations
if opt.const == 1 % Add constant to matrix of regressors
    XX = [XX ones(size(XX,1),1)];
end
% Add exogenous variables to matrix of regressors
XX = [XX exog(opt.p+1:end,:)]; 

n = size(YY,2); % Number of variables
ni = length(opt.ivar); % Number of variables of interest
nj = length(opt.jshock); % Number of shocks of interest
m = size(XX,2); % Number of parameters in each equation for y_t
nExog = opt.const + size(exog,2); % Number of exogenous variables
T = length(YY); % Number of observations used in estimating VAR

% Adjust indices representing narrative restrictions to account for
% losing the first p observations. If the restrictions are imposed in the 
% first p periods of the original dataset, drop the restrictions.
if ~isempty(restr.shockSignRestr)
    restr.shockSignRestr(:,3) = restr.shockSignRestr(:,3) - opt.p;
    restr.shockSignRestr(restr.shockSignRestr(:,3) <= 0) = [];
end
if ~isempty(restr.hdSignRestr)
    restr.hdSignRestr(:,3) = restr.hdSignRestr(:,3) - opt.p;
    restr.hdSignRestr(restr.hdSignRestr(:,3) <= 0) = [];
end

%% Conduct posterior inference on impulse responses.
% j-th column of B gives MLE of coefficients in j-th reduced form VAR 
% equation for y_t. Last rows correspond to exogenous terms (if included).
phiHat.B = (XX'*XX)\XX'*YY;
phiHat.S = (YY - XX*phiHat.B)'*(YY-XX*phiHat.B);
phiHat.Sigma = (1/T)*phiHat.S; % MLE of innovation covariance matrix
phiHat.P = (XX'*XX)\eye(m);
phiHat.cholP = chol(phiHat.P,'lower');

% Storage arrays.
rMinPost = zeros(opt.draws,opt.H+1,ni,nj);
rMaxPost = zeros(opt.draws,opt.H+1,ni,nj);
rSinglePriorPost = zeros(opt.draws,opt.H+1,ni,nj);
B = zeros(n*m,opt.draws);
Sigma = zeros(n,n,opt.draws);

phiDraw = 0; % Counter for no. of draws from posterior of phi
nonStableCount = 0; % Counter for no. of draws with nonstable VAR
draw = 0; % Counter for no. of draws with non-empty IS
Qempty = [];

while draw < opt.draws

    % Posterior sampler with independent improper (Jeffreys) prior
    phi.Sigma = iwishrnd(phiHat.S,T-m);
    phi.Sigmatr = chol(phi.Sigma,'lower');
    phi.Sigmatrinv = phi.Sigmatr\eye(n);
    phi.B = phiHat.B(:) + kron(phi.Sigmatr,phiHat.cholP)*randn(m*n,1);

    % Generate coefficients in orthogonal reduced-form VMA representation 
    % and check stability of VAR.
    [phi.vma,nonStable] = genVMA(phi,opt,nExog);
    nonStableCount = nonStableCount + nonStable;
    
    % If VAR is stable, proceed. Otherwise, redraw phi.
    
    if nonStable == 0
        
    phiDraw = phiDraw + 1;
    
    % Compute reduced-form VAR innovations.
    restr.U = (YY - XX*reshape(phi.B,m,n))';
        
    % Use Algorithm 1 (Step 2) to determine whether Q(phi|N,S) is nonempty
    % and draw Q satisfying restrictions.
    [Q0,restr.sphi,Qempty(phiDraw)] = drawQ0_ar(restr,phi,opt);
   
    if Qempty(phiDraw) == 0 % If Q(phi|N,S) is non-empty
            
            % Store reduced-form parameters.
            draw = draw + 1;
            Sigma(:,:,draw) = phi.Sigma;
            B(:,draw) = phi.B;

            % Compute draw from posterior of impulse responses given 
            % conditionally uniform prior for Q|phi.
            
            for hh = 1:opt.H+1 % For each horizon
    
                % Compute impulse response.
                rSinglePriorPost(draw,hh,:,:) = ...
                    phi.vma(opt.ivar,:,hh)*Q0(:,opt.jshock);

            end
            
           % Find upper and lower bounds of conditional identified set for
           % impulse responses.
           [rMinPost(draw,:,:,:),rMaxPost(draw,:,:,:),etaDraw] = ...
               approximateBounds(restr,phi,opt);

           clear etaDraw;

           if mod(opt.draws-draw,opt.dispIter) == 0
              
               fprintf('\n%d draws with non-empty identified set remaining...',...
               opt.draws-draw);
               
           end

    end
    
    end
    
end

% Compute posterior mean under single prior.
rSinglePriorMean = permute(mean(rSinglePriorPost,1),[2 3 4 1]);

% Compute lower and upper bound of set of posterior means.
meanlb = permute(mean(rMinPost,1),[2 3 4 1]);
meanub = permute(mean(rMaxPost,1),[2 3 4 1]);

% Compute robustified credible regions.
[credlb,credub] = credibleRegion(rMinPost,rMaxPost,opt);

% Compute highest posterior density (HPD) interval under single prior.
[hpdlb,hpdub] = highestPosteriorDensity(rSinglePriorPost,opt);

postMeanBoundWidth = meanub - meanlb; % Width of posterior mean bounds.
hpdWidth = hpdub - hpdlb; % Width of highest posterior density regions
credWidth = credub - credlb; % Width of robustified credible region

% Rows in Table 4 and 5 for this model:
disp('Table 4')
disp([meanlb(1,2) meanub(1,2) postMeanBoundWidth(1,2) rSinglePriorMean(1,2) max(abs(rSinglePriorMean(1,2)-meanlb(1,2)),abs(meanub(1,2)-rSinglePriorMean(1,2)))])

disp('Table 5')
disp([credlb(1,2) credub(1,2) credWidth(1,2) hpdlb(1,2) hpdub(1,2) hpdWidth(1,2) max(abs(credlb(1,2)-hpdlb(1,2)),abs(credub(1,2)-hpdub(1,2)))])
